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Abstract 

We describe a finite inhomogeneous three dimensional system of classical particles 
which interact through short and (or) long range interactions by means of a simple 
analytic spin model. The thermodynamic properties of the system are worked out 
in the framework of the grand canonical ensemble. It is shown that the system 
experiences a phase transition at fixed average density in the thermodynamic limit. 
The phase diagram and the caloric curve are constructed and compared with nu- 
merical simulations. The implications of our results concerning the caloric curve 
are discussed in connection with the interpretation of corresponding experimental 
data. 
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1 Introduction 



Up to now, our knowledge about the properties of highly excited nuclei generated by 
means of energetic particle or heavy ion collisions is still rather limited. In practice, the 
only well-established information concerns universal scaling properties [Q of fragment size 
multiplicities which are observed in the final stage of the reaction. These observables can 
be reproduced by means of simple, classical microscopic models. They can be interpreted 
as a sign for the existence of a phase transition of percolation type. 

One major handicap towards further understanding of excited nuclear matter is con- 
cerned with the distribution of the available energy among the nucleons in the excited 
system. At present it is not possible to say whether the highly excited expanding system 
can be considered to be in some thermal equilibrium or not. We know from specific exper- 
imental results that the nucleus is able to pick up and distribute huge amounts of energy, 
more or less equally, over internal degrees of freedom, at least in specific experimental 
circumstances. Theoretical interpretations in terms of so called dynamic models show 
that this is however not necessarily fully established. 

If one nevertheless assumes that equilibrium is reached, the central point of interest 
concerns the thermodynamic properties of finite and infinite nuclear matter up to very 
high energies and its connection with the universality properties of fragment distributions. 
During the last years several interesting attempts have been performed in this direction 
on the theoretical side. They rest on generic models, the Ising |J, the lattice gas and 
related models [4-7] and mean field approaches ||. They constitute interesting attempts 
towards a unified description of the global thermodynamics and fragment generation in 
nuclear systems. 

The work which we develop here takes place in the same framework. In practice we 
work out a classical, fully analytic realistic 3-dimensional model in the framework of the 
grand canonical ensemble. The analytic treatment allows for a clear insight into the 
correlations between the microscopic structure and the global thermodynamic properties 
of finite systems and their thermodynamic limit. It is able to incorporate both short and 
(or) long range interactions. Hence it can be used to take the Coulomb interaction into 
account and to investigate the thermodynamic properties of other classical microscopic 
systems. Finally many-body correlation functions which allow to define fragment sizes 
can be worked out analytically. 

In the sequel we shall essentially develop the expression of the partition function and 
its properties, and work out the relations between energy, temperature and density of the 
system. Indications about the determination of correlation functions will be given at the 
end of the letter and developed in the future. 

2 The model: grand canonical partition function for 
a system with short range interactions 

We consider a 3-dimensional system of classical particles confined in a volume V = (Nd) 3 
where N 3 is the number of cubic cells of linear dimension d (the extension to a parallelepi- 
pedic volume V = NiN 2 N 3 d 3 is straightforward). Each cell k is occupied by a particle 
(sk = 1) or empty (s^ = 0). Nearest neighbour particles interact by means of an attractive 
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2-body potential of strength — Vq (Vq > 0). The Hamiltonian of the system reads 



k=l <n-n>i 

where A is the total number of particles and < n ■ n- >i stands for "nearest-neighbour 
sites in each space direction. 

We shall show that this is the inhomogeneous counterpart of the homogeneous analytic 
model which we developped and discussed in ref. [ 20| . The model (1) can be brought into 
a form which looks more closely like the Ising model by introducing the spin variable 

a k = 2s k -l (2) 

Up to the kinetic term, Hq then becomes an Ising Hamiltonian in a magnetic field 

A N 3 

#o = $>fc/2m (V /4)[ J2 cT k a £ + 6j2^ + SN 3 ] (3) 

k=l <n-n>i k=l 

The conservation of the number of particles A is expressed through the constraint 

N 3 

(1/2) £(1 + a k ) = A (4) 

k=i 

In order to get an exactly solvable model we replace the Ising-like model (3) by the 
spherical model ||. It has been known since a long time that the spherical model is 
exactly solvable in any number of dimensions, in the presence of symmetry-breaking 
fields, for very general (including long-range) interactions and even in the presence of 
competing interactions or disorder [9-17]. 

Yet it does display, in more than two and less than four dimensions, a second order 
phase transition with non-trivial (i.e. non mean field) critical exponents. The spherical 
model represents therefore an ideal test bench on which complicated concepts can be 
checked analytically. 

To introduce it, we replace the Ising spins a k = ±1 in (3) by continuous real variables 
o k G [— oo, +oo] and constrain these variables through the relation 

N 3 

E^ = Ar3 ( 5 ) 

k=l 

Practically it is easier to implement the constraints (4), (5) only in the mean through 
Lagrange multipliers A and fi. In the thermodynamic limit iV — > oo it is known that this 
leads to the same result JTBJ. The mean spherical Hamiltonian can then be written as 

N 3 N 3 

H = H + \J24 + (a*/2)(XI °* + N ") ( 6 ) 

k=l k=l 

and the partition function becomes 

/+oo f+oo 
Y[dpl l[da k exp(-{3H). (7) 

-°° kA J -°° k 
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We remark that going over from the Ising model to the spherical model implies the 
replacement 



k=l{a k =±l} k=l 



+00 



N 3 N 3 

n e - n / ''"<~ 



00 



The quantity p\ is the momentum of particle k in the space direction i (i = 1, . . . , n). In 
the sequel we fix n = 3 and assume periodic boundary conditions in the 3 space directions. 
The constraints 

-P^dlnZ/dX = N 3 

-(3- 1 din Z/dfi = A (8) 

fix the parameters A and \x. This completes the definition of the model. 

At first sight the spherical constraint may appear bizarre, introducing an effective long- 
range interaction between the spins a^. But this is only apparent and can be interpreted 
as follows. Consider an m-component spin variable S a , a = 1, ... ,m of unit length 

m m 

S a S a = 1. The 0{m) vector model with Hamiltonian H ( m ) = —J Yl Yl ^i^j can 

a=l <n n>i a=l 

be shown to converge to the spherical model in the limit m —>■ 00 ]10[ (for m = 1,2,3 one 
recovers the familiar Ising, XY and Heisenberg models respectively). 

In our case, the procedure of replacing the model (1) by the spherical model (6-8) 
leads to an extension which describes the fact that the particles do not lie at constant 
distance from each other, since their centres in cells of volume d 3 can be found at random 
inside each occupied cell. This is in fact the essence of an earlier modelisation of the 
system pOf . Allowing Sk (resp. a^) to vary continuously can be interpreted in practice 
as the description of a system for which = ±1 and an interaction V <n . n>1 which 
takes values on the interval [0, +00]. But events with large value of V <n . n>1 are strongly 
suppressed because the largest part of their contributions are generated by a' k s lying in a 
strip of values, due to the quadratic form of H in (7) and the constraint (5) on the volume. 
Our interpretation does not take care of the fact that the positions of a set of interacting 
particles are correlated and hence, so are the two-body strengths V <n . n>1 . However, we 
expect that this is not a major problem because of the fact that our calculations concern 
the partition function Z which delivers average quantities and hence these correlations 
should not show a strong effect. 

In order to diagonalise the free energy in the exponent of Z we Fourier-transform {<Jk\ 
to {tj}. In each space direction 

N 

T£ = JM 2^ fc ' 



= 1 



Then 



k I 
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The quadratic form in {r^} can be explicitly integrated and 

Z(J3, A, N 3 ) oc (2vrm) 3A/2 • /r^^/ 2 • exp(-/3/iiV 3 /2 + 3(3V N 3 /A) 

■ exp {(3 [(// - 3V )/2] 2 N 3 / [4(A - 3V /2)}} 

■ Y[ [A-K)(coseu+cos^ 2 + cos^ 3 )/2]- 1/2 (10) 

{*1,*2,«3} 

with = 2nik/N (i^ — 0, 1, . . . , N — 1) and the product runs over all possible values 
of (ii,*2 5 *3)- The first term on the r.h.s. of (10) corresponds to the kinetic energy 
contribution. 
We define 

M = -[1 + 2((3N 3 )- 1 d\nZ/d/2} = -(fi - 3V )/[A(X - 3V /2)} 
= -1 + 2A/N 3 

which is constant and related to the density through 

M = 2p/p - 1 

where po is the normal density of nuclear matter. Then equations (8) and (10) lead to 

N*-l 

]T [A(/3) - = 2iV 3 /3(l - M 2 ) 

i=0 



with aa = a hi2 i 3 = cos £ ik 

k=l 

and n(P)=4{3V /2-\{/3)]M + 3Vo. (11) 

Equation (11) is a kind of dispersion relation with N 3 poles. We look for the solution 
which corresponds to A > Voa; max /2 which is positive (a max = 3) and guarantees that the 
integrals in the expression of Z make sense. 

From (11) one can immediately see that in the thermodynamic limit (N — > oo) 
A(/3) — > V a max /2. This shows a singularity in the total free energy, even though the free 
energy per particle remains finite. 

For M = (p/po = 1/2), this transition is continuous and occurs at the critical 
temperature Jl3[ 

(3 C = 2W 3 (0)/V 

POO 

where W 3 (0) = (1/2) / e - 3x I 3 (x)dx = 0.2527 

Jo 



Iq(x) being a modified Bessel function |2T 



It is now possible to construct the caloric curve for fixed density. The energy per 
particle is given by 



3 V d\nZ 3 Jr f 3N 3 r „ „ „ „ „, 

-— ^ — = V Q < \2M(M + 1) + 1 

2/5 A 09 2(3 \ AA 1 K ' J 



iV 3 -l 



+ 4A S (A/3 -0^/2)1 (12) 



with 9 = /3Vq. The first term in the expression corresponds to the kinetic energy contri- 
bution. 

Contrary to the model introduced in ref [|13j one sees that the solutions (A, jj) obtained 



through relations (11) are of the same type for any value of M because here M is a constant 
quantity. Hence one concludes that the phase transition is of the same order for any fixed 
density. 

The energy e given by (12) is a continuous function of (3. The specific heat Cy reads 



C v = -(3 2 (de/d(3) = 3/2 -(V Q /4A)[j2 

i=0 



N 3 -l 



(A - Vooti/2) 



This quantity is again continuous in (3. It is in fact well known []2"2"|, 01 that the system 



shows a phase transition at (3 C (M) where Cy ~ ((3 C — (3)~ a and a = — 1. There appears 
however a discontinuity in the derivative of Cy with respect to (3 in the thermodynamic 
limit. This can indeed be guessed numerically in Fig. 1 for increasing sizes of the system. 
We also remark that the convergence towards the N — > oo limit appears to be relatively 
rapid, even around the critical point. For (3 > /3 C (M) the specific heat is constant. 



3 Extension of the description to a system with short 
and (or) long range interactions 

The model can be easily extended to the case where particles interact by means of short 
and (or) long range two-body interactions. Hence one can for instance investigate the 
effects of the Coulomb interaction in the excited nuclear system. 
The Hamiltonian can be written in the general form 

H = K + V = K-[V 01 SiSj + Vw s i s J + --- 

<n-n.>i <n-n.>2 

■■■ + v 0n s * s j}+ ^J2 Si + x J2( 2si ~ 1 ^ 2 ( 14 ) 

<n-n> n i i 

where K is the kinetic energy and < n-n. >k stands for k th nearest neighbour in all space 
directions and 

V 01 = V s + V e {1) 
V 02 = V} 2) 



V — T/( n ) 
VOn — V i 

with Vs > and vf 1 < (i = 1, . . . , n). Vs stands for the short range potential and V^ 
for the long (Coulomb) contribution. We parametrise this strength in the following way 

V e {k) = (Z p /A)e 2 /kd (15) 
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where Z p is the total number of protons. Here again we introduce periodic boundary 
conditions on the variables <7j = 2s« — 1. 

The present procedure smears out the charge over all existing particles. Hence the 
interaction is not taken into account exactly, but in the average over the occupied parts 
of the volume, since protons are not explicitly distinguished from neutrons. This should 
not affect strongly the calculations since one expects that the charges are effectively dis- 
tributed approximately uniformly in the ratio of Z p /A over the set of occupied cells. It 
should also be mentioned that, as in the case of the short range interaction, the particles 
interact only along the three space directions and not the directions of the diagonals. 
These contributions can be approximately taken into account by a consequent renormal- 
ization of the interaction strength. 

The grand canonical partition function can now be worked out by following the same 
procedure as before. Notice that the extension of the spin variables {<Xfc} to the continuum 
limit introduces, as previously in section 2, values of |Vbfc|, in the interval [0,+oo]. If 
one follows the interpretation developed there, the long range interaction in (13) can in 
principle take large values. However, like there, large values of Vok are suppressed because 
of the quadratic form of H and the constraint (5). After some algebra one gets 

Z(0, A, N 3 ) oc {2<nm) 3A l 2 • /T (3A+iv3)/2 • exp(-/^7V 3 /2 + 3f3W N 3 /A) 

■ exp {(3 [(/i - 3W )/2} 2 N 3 /4(X - 3W /2) } 

-1/2 



n 

{n,«2,«3} 



A - J2 K^cos^ + cos#> + cos^)/2 

p=l 



(16) 



with = %pi e (i e = 0,l,...,N-l),(p=l,...,n) and W = £ V 0p 

P =i 

In the same way as before, A(/3) and are obtained through 

= 2N 3 f3(l - M 2 ) 



E 

{n,«2,«3} 



HP) - E ^(cosdf + cosd 2 P) + cosC ( 3 P) )/2 
P =i 

M = -1 + 2A/N 3 

and fi((3) = 4(3W /2 - \{(3))M + 3W (17) 



The energy per particle e(f3) is the same as the one given by (12) with Vq replaced by Wo, 

3 X ^ c) In Z 

e =2-p-A^ Vop ^t With B * = PV <* 
p=i y 

and an expression similar to (13) for the specific heat Cy. 

4 Results 

The typical behaviour of the caloric curve for different densities is shown in Figs. 2-4 
for the case of an attractive short range, an attractive long range and a combination of 
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an attractive short range with a repulsive long range interaction respectively. In each 
case, the energy increases with respect to the temperature and one does not observe the 
presence of a plateau, even though the description predicts a phase transition at every 
density. One observes simply some more or less pronounced kink in the slope of the curve 
which corresponds in practice to the temperature at which the transition occurs. This 
is observed for all types of 2-body interactions. As already mentioned in section 3, the 
reason for this smooth behaviour is related to the fact that in the spherical model the 
specific heat remains finite at the transition point even at the thermodynamic limit. Only 
the derivative of the specific heat shows a discontinuity at the transition points. From the 
knowledge of these transition points it is possible to construct the phase diagram in the 
(T/T c , p/po) plane, Fig. 5. Here T c corresponds to the critical point at p/po = 0.5 where 
the usual spherical model undergoes a second order transition working at zero magnetic 
field h ||13|| . In the spherical model, the transition for p/po ^ 0.5 is of first order, owing 
to the discontinuity in the curve h vs M at T < T c . This means that when h = 0, the 
system can pass from a magnetization +M to — M. But in the present model we fix M to 
a certain value (constraint (8)), so flips from +M to — M are not possible here. Therefore, 
one expects a behaviour for every M similar to what happens at M = 0. This is why we 
obtain the same smooth behaviour at every value of p/po- 



5 Numerical simulations in the framework of a cel- 
lular model 

The present description is close to former numerical approaches (cellular model) which 
were introduced in order to study the thermodynamic properties and cluster content of an 
excited disordered system of particles |]n| |2U|. However in these simulations the system 
was chosen to be homogeneous (p/po = 1) an d the numerical procedure did not imply 
that the system was in thermodynamic equilibrium. There was no sign for the existence 
of a thermodynamic phase transition. 

We present here the results of numerical simulations which generalise in the framework 
of the canonical ensemble the former calculations. We consider a finite volume of space 
V in 3D made of unit cells of linear dimensions d = 1.8 fm, V = N 3 .d. Each cell is either 
occupied or unoccupied, the total number of particles is A such that the normalized 
density p/po = A/N 3 . The particles interact by means of a nuclear [23[] and Coulomb 
interaction. 

For fixed density and temperature the equilibrium energy of the system is obtained 
by means of the following algorithm [J24|. Assign an initial configuration with A particles 



in N 3 cells. Select randomly with equal probability one of the two following procedures. 
In the first procedure one selects randomly a particle among the A occupied cells and 
chooses a new position of this particle inside its cell. This changes the energy of the 
system by an amount AE which can be positive or negative. At this stage one applies 
the Metropolis Monte Carlo algorithm before returning to the procedure selection step. 
In the second procedure one selects randomly a particle "a" in one of the occupied cells 
and an empty cell "6" among the N 3 — A empty cells. One creates a particle at the center 
of the cell "6" with the same charge (proton or neutron) as the charge of the particle 
at "a". One destroys the particle at "a" and hence guarantees total mass and charge 
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conservation. This changes the energy of the system by an amount AE and one applies 
again the Metropolis Monte Carlo algorithm before going back to the procedure selection 
step. The average number of measurements is 4000 corresponding to ten times more 
Metropolis steps. 

Typical results are shown in Fig. 6. The fluctuation of the results is not larger than the 
size of the dots. The energy per particle in the abscissa is the sum of the interaction and 
the kinetic energy which gives the trivial contribution 3T/2. We checked the consistency 
of the procedure by means of numerical simulations of the kinetic energy of a system of 
free particles from a Boltzmann distribution. The fluctuations in the kinetic energy are 
negligible. 

The behaviour of the caloric curve is very similar to the curve obtained by means of 
the analytic model, see Fig. 4. Indeed, the line which goes through the points corresponds 
to a calculation with the spherical model and Vq = 5.5 MeV for the attractive short range 
interaction. It is worthwile to notice that there is no sign for the existence of a phase 
transition. If however we interpret the results in the framework of the spherical model, the 
observed smooth behaviour would correspond indeed to a continuous transition reflected 
on the discontinuity in the derivative of the specific heat at the temperature predicted by 
the model. 



6 Discussion and final remarks 

The fact that the description leads to a continuous phase transition can a priori seem 
somewhat surprising. Indeed, other models which are related to spin systems like the 
lattice gas model clearly show a first order transition [4-7] which is not reproduced here. 
One could argue that we have 'softened' the strength of the transition by working in the 
framework of the spherical model, because it corresponds to the limit m — > oo in the 0(m) 
model. However, we remark that this is not the case. We obtain a continuous transition 
owing to the constraint fixing the magnetization, which prevents it from having the typical 
fluctuations which happen at the first order transition points in these other models. In 
fact, one can consider this constraint in the Ising model and see numerically that the whole 
transition line corresponds to a second order phase transition [B3|. We have considered 



here the spherical model approximation in order to obtain analytical expressions for the 
different quantities and to incorporate short and long range interactions. The controversy 
between the order of the transition should work as an incentive to experimental physics to 
look for the best way to identify the specific properties of excited nuclei and nuclear matter 
and try to confirm or infirm the existence of a discontinuous phase transition [26-29]. 

The numerical simulations on the inhomogeneous model have been made in the frame- 
work of the canonical ensemble, with strictly fixed volume and number of particles. The 
fact that the results are in qualitative agreement with those obtained analytically in the 
framework of the grand canonical ensemble indicates that volume and particle fluctua- 
tions do not strongly affect the physical issue. It also shows that the interpretation of the 
continuum limit approximation leading to the spherical model description makes sense. 
Notice that in the experiment, volume and exact particle number can hardly be rigorously 
fixed. It would nevertheless be interesting to work out the dynamics in the framework of 
the microcanonical ensemble where the isolated nuclear system is treated rigorously |3C 
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Whether this is possible has to be investigated in practice. 

Finally it would be of particular interest to work out averaged n-body correlation 
functions which allows to determine the fragment size distributions. This can in principle 
be done by means of a straightforward extension of the partition function Z to a generating 

function by adding terms of the type exp X^i(^i s « + <?i(l ~~ s i))j into the sum over 
occupation numbers. Here hi and gi are arbitrary real constants and the probability 
of an isolated cluster with (ii, . . . ,i n ) occupied cells can be obtained by deriving this 
generating function with respect to specific parameters hi and (which fix the surface 
of the cluster) and putting the parameters to zero. The determination of the cluster 
formation probabilities would then lead to a complete description of the fragmentation of 
the considered system in the (e, p, f3) space. 
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Figure caption 



Figure 1: Specific heat calculated from (13) with p/p = 0.90, V = 8MeV, and different 
sizes of the system (N = 10, 20, ... , 50). 

Figure 2: Caloric curve obtained with an attractive short range (nearest-neighbour) in- 
teraction with strength V = 8MeV, A = N 3 p/p = 250. 

Figure 3: Caloric curve obtained with an attractive long range interaction (n = 20) with 
strength W = 8 MeV. 

Figure 4: Caloric curve for systems which experience both an attractive short range 
(nearest-neighbour) with strength Vq = 8 MeV and a repulsive long range interaction 
with strengths given by (15) (n = 20). Here Z p /A = 0.394 and d = 1.8 fm, see eq.(14). 

Figure 5: Phase diagram T/T c vs. p/po- The transition line is read from the behaviour 
of the specific heat. The point T/T c = 1 corresponds to the continuous phase transition 
in the thermodynamic limit of the usual spherical model. 

Figure 6: Caloric curve obtained by means of Metropolis Monte Carlo simulations. The 
points correspond to numerical results for N 3 = 216, p/po = 0.5, A = 108 and total 
charge number Z p = 47. The full line corresponds to the results obtained through the 
spherical model with V = 5.5 MeV. The arrow labelled T indicates the phase transition 
temperature. See explanations and comments in the text. 
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Caloric curve for attractive short range interaction 



1 1 1 1 1 1 1 1 1 1 

p/p =0.30 / 


i 

p/po=^ 


p/p =0J0 / 

\ a i i i i i i i i i 


p/p =0.90 / 

r 1 1 1 1 1 1 1 1 1 



40 -20 20 40 -20 20 40 

Energy (MeV) 

Fig. 2 



Caloric curve for attractive long range interaction 
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Fig. 4 
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